!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Channel test case
!!
!! \n
!!
!! We assume the following boundary labels:
!! <ul>
!!  <li>left inflow: 3
!!  <li>right outflow: 2
!!  <li>upper boundary: 1
!!  <li>bottom boundary: 5, 4
!!  <li>bottom obstacle: 6, 7
!! </ul>
!<----------------------------------------------------------------------
module mod_channel_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_channel_test_constructor, &
   mod_channel_test_destructor,  &
   mod_channel_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof,         &
   coeff_dir,        &
   coeff_norm,       &
   coeff_sponge,     &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter   :: test_name = "channel"
 character(len=100), protected :: test_description(3)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 0
 character(len=*), parameter   :: ref_prof = "none"

 logical, protected :: mod_channel_test_initialized = .false.
 ! private members
 real(wp) :: &
   nu,     & ! viscosity
   t_i,    & ! inflow temperature
   p_i,    & ! inflow pressure
   u_i,    & ! inflow velocity
   ! obstacle geometry
   x_c, z_c, r_x, r_z
 character(len=*), parameter :: &
   test_input_file_name = 'channel.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_channel_test'

 ! Input namelist
 namelist /input/ &
   nu, u_i, x_c, z_c, r_x, r_z

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_channel_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_channel_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Channel, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,';'
   write(test_description(3),'(a,e9.3,a)') &
     '  input velocity = ',u_i,'.'

   ! This test case does not consider gravity nor Earth rotation
   phc%omega   = 0.0_wp
   phc%gravity = 0.0_wp

   ! Define the input state
   p_i = phc%p_s
   t_i = 300.0_wp

   mod_channel_test_initialized = .true.
 end subroutine mod_channel_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_channel_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_channel_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_channel_test_initialized = .false.
 end subroutine mod_channel_test_destructor

!-----------------------------------------------------------------------

 !> Uniform flow using the input values.
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

   u0(1,:) = p_i/(phc%rgas*t_i)                       ! rho
   u0(2,:) = u0(1,:) * ( phc%cv*t_i + 0.5_wp*u_i**2 ) ! e
   u0(3,:) = u0(1,:) * u_i                            ! Ux
   u0(4:,:) = 0.0_wp                                  ! Uyz
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

 end function coeff_init
 
!-----------------------------------------------------------------------

 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1),size(x,2))

   u0 = 0.0_wp

 end function coeff_outref

!-----------------------------------------------------------------------

 !> Use the initial condition.
 pure function coeff_dir(x,u_r,breg) result(ub)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: ub(2+size(x,1),size(x,2))

   ub = coeff_init(x,u_r)

 end function coeff_dir

!-----------------------------------------------------------------------

 !> Normal to the boundary
 !!
 !! \warning This function should be only called for the bottom
 !! boundary.
 !!
 !! In 2D, we have
 !! \f{equation}{
 !!  y(x) = z_c + r_z\sqrt{1-\left(\frac{x-x_c}{r_x}\right)^2}
 !! \f}
 !! so that
 !! \f{equation}{
 !!  \tilde{\underline{n}}(x) = 
 !!   \left[\begin{array}{c} \frac{dy}{dx}(x) \\\ -1 \end{array} \right]
 !!   = \left[\begin{array}{c} 
 !!    -\frac{r_z}{r_x}\frac{\xi}{\sqrt{1-\xi^2}} \\\ -1
 !!   \end{array} \right]
 !! \f}
 !! where the \f$\tilde{}\f$ indicates that the resulting vector must be
 !! normalized and
 !! \f{equation}{
 !!  \xi = \frac{x-x_c}{r_x}.
 !! \f}
 !<
 pure function coeff_norm(x,breg) result(nb)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: nb(size(x,1),size(x,2))

  integer :: l
  real(wp) :: xi(size(x,2))

   breg_case: select case(breg)
    case(4,5) ! bottom boundry, 2D
     nb(1,:) =  0.0_wp
     nb(2,:) = -1.0_wp
    case(6,7) ! obstacle in the bottom boundary
     xi = (x(1,:) - x_c)/r_x
     nb(1,:) = -r_z/r_x * xi/sqrt(1.0_wp-xi**2)
     nb(2,:) = -1.0
     do l=1,size(nb,2)
       nb(:,l) = nb(:,l)/sqrt(nb(1,l)**2+nb(2,l)**2)
     enddo
    case default
     nb = huge(1.0_wp) ! this is an error
   end select breg_case

 end function coeff_norm

!-----------------------------------------------------------------------

 pure subroutine coeff_sponge(tmask,tau,taup,x)
  real(wp), intent(in) :: x(:,:)
  logical, intent(out) :: tmask
  real(wp), intent(out) :: tau(:), taup(:)

  real(wp), parameter :: x_r =  2.5_wp
  real(wp), parameter :: dx = 0.5_wp
  real(wp), parameter :: big = huge(1.0_wp)
  integer :: l
  real(wp) :: d_bnd, coeff

   ! initialization
   tau = 0.0_wp
   taup = 0.0_wp
   tmask = .false.

   do l=1,size(x,2)
     d_bnd = (x_r-x(1,l))/dx
     if(d_bnd.le.3.0_wp) then
       tmask = .true.
       if(d_bnd.lt.1000.0_wp/big) then
         coeff = big/1000.0_wp
       else
         coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
       endif
       tau(l)  = coeff
       taup(l) = coeff
     endif
   enddo

  end subroutine coeff_sponge

!-----------------------------------------------------------------------

end module mod_channel_test

